#!/usr/bin/perl

#TODO:
#1. if amount of latency in group is close, don't worry about off by 1 or 2 buckets

sub usage {
	print "usage: stat_analysis.pl -f1 <file> -f2 <file> -p1 -p2 -c <command> -debug -h\n";
	print "-f1: file containing first set of profiles\n";
	print "-f2: file containing second set of profiles\n";
	print "-p1: do NOT do first phase of filtering\n";
	print "-p2: do second phase of filtering\n";
	print "-c: command to run on the data.\n";
	print "    acceptable commands:\n";
	print "\tchisquare : get the chi-square value\n";
	print "\tgrouplat: compare group positions and latencies (must also use -p2 option)\n";
	print "\tgroupops: compare group positions and operations (must also use -p2 option)\n";
	print "\ttotlat: compare total latencies\n";
	print "\ttotops: compare total operations\n";
	print "\tcorrelation: statistical correlation\n";
	print "\tearthmover: earth mover's distance (EMD)\n";
	print "\tquadratic: quadratic-form distance\n";
	print "\tnone: just do phase 1 and/or 2\n";
	print "-debug: print debug information\n";
	print "-h: print this help menu and exit\n";
	exit 0;
}

sub max {
	my $m = shift @_;
	foreach (@_) {
		if ($_ > $m) {
			$m = $_;
		}
	}
	return $m;
}

sub min {
	my $m = shift @_;
	foreach (@_) {
		if ($_ < $m) {
			$m = $_;
		}
	}
	return $m;
}

sub get_totlat {
	my @data;
	my $tot = 0;
	my $fname = $_[0];
	open INFILE, "<$fname" or die "Cannot open file: $!";

	while (<INFILE>) {
		chomp();
	
		if (/(^OP_[\w]+) /) {
			@data = split / /, $_;
			$tot += $data[2];
		}
	}

	close INFILE;
	return $tot;
}

sub groupanal {
	my $i, $numgroups, $start, $end, $peak, $increase, $grouptot, $bucketweight, $groupweight, @res;

	$start = 0;
	$numgroups = 0;
	$bucketweight = 32;
	for ($i = 3; $i <= $#_; $i++, $bucketweight*=2) {
		if (!$start) {
			if ($_[$i]) {
				$start = $i - 2;
				$peak = $_[$i];
				$increase = 1;
				$grouptot = $_[$i];
				$groupweight = $_[$i] * $bucketweight * (3/4);
			}
		} else {
			if (!$_[$i]) {
				$end = $i - 3;
				if ((($_[$i - 1] < 10) && ($_[$i - 1] > 0)) && (($_[$i + 1] < 10) && ($_[$i + 1] > 0))) {
					
				} elsif ((($grouptot < ($_[1]/1000)) && ($groupweight < ($_[2]/20))) || (($groupweight < ($_[2]/1000)) && ($grouptot < ($_[1]/20)))) {
					if ($debug) {
						print "ignoring group: $start - $end (peak: $peak, tot: $grouptot, weight: $groupweight)\n";
					}
					$start = 0;
				} else {
					$numgroups++;
					$res[0] = $numgroups;
					$res[$numgroups] = {};
					$res[$numgroups]->{START} = $start;
					$res[$numgroups]->{END} = $end;
					$res[$numgroups]->{GROUPWEIGHT} = $groupweight;
					$res[$numgroups]->{GROUPTOT} = $grouptot;
					if ($command eq "grouplat") {
						$res[$numgroups]->{AMOUNT} = $groupweight;
					} else {
						$res[$numgroups]->{AMOUNT} = $grouptot;
					}
					if ($debug) {
						print "group: $start - $end (peak: $peak, tot: $grouptot, weight: $groupweight)\n";
					}
					$start = 0;
				}
			} elsif ($increase) {
				if ($_[$i] > $peak) {
					$peak = $_[$i];
				} elsif ((length($_[$i]) > 1) && (length($_[$i]) < length($peak))) {
					$increase = 0;
				}
				$grouptot += $_[$i];
				$groupweight += $_[$i] * $bucketweight * (3/4);
			} else {
				if ((length($_[$i]) > length($_[$i - 1])) && ((length($_[$i - 1]) > 1) || (length($_[$i] > 2)))) {
					$increase = 1;
					$end = $i - 3;
					$numgroups++;
					$res[0] = $numgroups;
					$res[$numgroups] = {};
					$res[$numgroups]->{START} = $start;
					$res[$numgroups]->{END} = $end;
					$res[$numgroups]->{GROUPWEIGHT} = $groupweight;
					$res[$numgroups]->{GROUPTOT} = $grouptot;
					if ($command eq "grouplat") {
						$res[$numgroups]->{AMOUNT} = $groupweight;
					} else {
						$res[$numgroups]->{AMOUNT} = $grouptot;
					}
					if ($debug) {
						print "group: $start - $end (peak: $peak, tot: $grouptot, weight: $groupweight)\n";
					}
					$start = $i - 2;
					$peak = $_[$i];
					$grouptot = $_[$i];
					$groupweight = $_[$i] * $bucketweight * (3/4);
				} else {
					$grouptot += $_[$i];
					$groupweight += $_[$i] * $bucketweight * (3/4);
				}
			}
		}
	}

	if ($start) {
		$end = $i - 3;
		$numgroups++;
		$res[0] = $numgroups;
		$res[$numgroups]->{START} = $start;
		$res[$numgroups]->{END} = $end;
		$res[$numgroups]->{GROUPWEIGHT} = $groupweight;
		$res[$numgroups]->{GROUPTOT} = $grouptot;
		if ($command eq "grouplat") {
			$res[$numgroups]->{AMOUNT} = $groupweight;
		} else {
			$res[$numgroups]->{AMOUNT} = $grouptot;
		}
		if ($debug) {
			print "group: $start - $end (peak: $peak, tot: $grouptot, weight: $groupweight)\n";
		}
	}

	return @res;
}

use Statistics::Distributions;

# Calculate the chi-square for two sets of bins.  A small value indicates a
# significant difference between the distributions.
sub chstwo {
	my @bins1, @bins2;
	my $j, $temp, $df, $chsq, $nbins, $highest;

	@bins1 = @{$_[0]};
	@bins2 = @{$_[1]};

	$nbins = $#bins1 - 2;
	$df = $nbins;
	$chsq = 0.0;

	for ($j = 3; $j <= $#bins1; $j++) {
		if (($bins1[$j] == 0.0) && ($bins2[$j] == 0.0)) {
			--$df;
		} else {
			$temp = length($bins1[$j]) - length($bins2[$j]);
			$chsq += ($temp * $temp) / (length($bins1[$j]) + length($bins2[$j]));
		}
	}	

	if ($debug) {
		printf("df: %d, chsq: %f\n", int(0.5 * $df), (0.5 * $chsq));
	}
	if (!int(0.5 * $df)) {
		return 0.0;
	} else {
		return &Statistics::Distributions::chisqrprob(int(0.5 * $df), (0.5 * $chsq));
	}
}

sub correlate {
	my @bins1, @bins2;
	my $j, $xmean, $ymean, $Sxx, $Syy, $xi, $yi, $xdiff, $ydiff, $sum;

	@bins1 = @{$_[0]};
	@bins2 = @{$_[1]};

	$xmean = $ymean = $Sxx = $Syy = $sum = 0;
	for ($j = 3; $j <= $#bins1; $j++) {
		if (!$bins1[$j] && !$bins2[$j]) {
			next;
		}
		$xmean += length($bins1[$j]);
		$ymean += length($bins2[$j]);
	}

	for ($j = 3; $j <= $#bins1; $j++) {
		if (!$bins1[$j] && !$bins2[$j]) {
			next;
		}
		$xi = length($bins1[$j]);
		$yi = length($bins2[$j]);
		$xdiff = $xi - $xmean;
		$ydiff = $yi - $ymean;
		$sum += ($xdiff * $ydiff);
		$Sxx += ($xdiff * $xdiff);
		$Syy += ($ydiff * $ydiff);
	}

	return ($sum/sqrt($Sxx*$Syy)) * 100;
}

sub earthmover {
	my @bins1, @bins2;
	my $tot1, $tot2, $j, $carrying, $emd, $max;

	@bins1 = @{$_[0]};
	@bins2 = @{$_[1]};

	$tot1 = $tot2 = 0;
	for ($j = 3; $j <= $#bins1; $j++) {
		$bins1[$j] = $bins1[$j] ? length($bins1[$j]) : 0;
		$bins2[$j] = $bins2[$j] ? length($bins2[$j]) : 0;
		$tot1 += $bins1[$j];
		$tot2 += $bins2[$j];
	}

	$carrying = $emd = 0;
	for ($j = 3; $j <= $#bins1; $j++) {
		my $a = $bins1[$j]/$tot1;
		my $b = $bins2[$j]/$tot2;
		$carrying += $a - $b;
		$emd += abs($carrying);
	}

	$max = ($#bins1 - 4);
	return ($emd/$max) * 100 * 5;
}

use Math::MatrixReal;

# Quadratic-form distance: d = sqrt(trans(h - k) * A * (h - k)) where h and k
# are vectors that list all the entries in H and K. Cross-bin information is
# incorporated via a similarity matrix A where aij denote
# similarity between bins i and j . Here i and j are sequential (scalar)
# indices into the bins.  For our experiments, we followed the recommendation
# in Niblack et al. (1993) and used aij = 1 - (dij/dmax) where dij is the
# ground distance between bins i and j of the histogram, and dmax = max(dij).

sub quadratic {
	my @bins1, @bins2;
	my $histsize, $i, $j, $tot1, $tot2, $dmax;
	my $m1, $m2, $A, $ans;

	@bins1 = @{$_[0]};
	@bins2 = @{$_[1]};

	$histsize = &max($#bins1, $#bins2);
	$m1 = new Math::MatrixReal(1, $histsize);
	$A = new Math::MatrixReal($histsize, $histsize);
	$m2 = new Math::MatrixReal($histsize, 1);
	$ans = new Math::MatrixReal(1, 1);

	$tot1 = $tot2 = 0;
	for ($j = 3; $j <= $histsize; $j++) {
		$tot1 += $bins1[$j];
		$tot2 += $bins2[$j];
	}

	for ($j = 3; $j <= $histsize; $j++) {
		$bins1[$j] = $bins1[$j]/$tot1;
		$bins2[$j] = $bins2[$j]/$tot2;
	}

	for ($j = 3; $j <= $#bins1; $j++) {
		$m1->assign(1, $j, $bins1[$j]);
	}

	for ($j = 3; $j <= $#bins2; $j++) {
		$m2->assign($j, 1, $bins1[$j]);
	}

	$dmax = 0;
	for ($i = 1; $i <= $histsize; $i++) {
		for ($j = 1; $j <= $histsize; $j++) {
			$A->assign($i, $j, abs($m1->element(1, $j) - $m2->element($i, 1)));
			$dmax = ($A->element($i, $j) > $dmax) ? $A->element($i, $j) : $dmax;
		}
	}

	for ($i = 1; $i <= $histsize; $i++) {
		for ($j = 1; $j <= $histsize; $j++) {
			$A->assign($i, $j, 1 - ($A->element($i, $j) / $dmax));
		}
	}

	$ans = $m1 * $A * $m2;
	return ((1 - sqrt($ans->element(1,1))) * 100);
}

##############################################################################
# main
##############################################################################

use Switch;

$phase1 = 1;
$phase2 = 0;

for ($argnum = 0; $argnum <= $#ARGV; $argnum++) { 
	switch ($ARGV[$argnum]) {
		case "-f1"	{ $filename1=$ARGV[++$argnum]; $op1 = 1; }
		case "-f2"	{ $filename2=$ARGV[++$argnum]; $op2 = 1; }
		case "-c"	{ $command=$ARGV[++$argnum]; $op3 = 1; }
		case "-p1"	{ $phase1 = 0; }
		case "-p2"	{ $phase2 = 1; }
		case "-debug"	{ $debug = 1; }
		case "-h"	{ &usage; }
		else		{ print "unknown arg: $ARGV[$argnum]\n"; &usage; }
	}
}

if (!$op1 || !$op2 || !$op3) {
	&usage;
}

if ($debug) {
	print "file1: $filename1, file2: $filename2, command: $command\n";
}

$totlat_file1 = &get_totlat($filename1);
$totlat_file2 = &get_totlat($filename2);

if (($command eq "grouplat") || ($command eq "groupops")) {
	$phase2 = 1;
}

open INFILE1, "<$filename1" or die "Cannot open file: $!";

IN: while (<INFILE1>) {
	chomp();
	if (/(^OP_[\w]+) /) {
		$op = $1;
		@data1 = split / /, $_;
	} else {
		next;
	}

	$found = 0;
	open INFILE2, "<$filename2" or die "Cannot open file: $!";
	while (<INFILE2>) {
		chomp();
		if (/(^$op) /) {
			@data2 = split / /, $_;
			$found = 1;
			last;
		} else {
			next;
		}
	}
	close INFILE2;

	if ($found == 0) {
		print "WARNING: Did not find $op in $filename2\n";
		next;
	}

	if ($debug) {
		print "\ndata1: @data1\ndata2: @data2\n";
	}

	# PHASE 1: sift out anything insignificant or very similar
	if ($phase1) {
		if ($debug) {
			print "\n";
		}

		if (($data1[2] < ($totlat_file1 / 100)) && ($data2[2] < ($totlat_file2 / 100))) {
			if ($debug) { print "latency is insignificant\n"; }
			print "$op: 0.0\n";
			next IN;
		}

		if ((($data1[2] < 5) && ($data2[2] < 5)) || (($data1[1] < 3) && ($data2[1] < 3))) {
			if ($debug) { print "number of operations is insignificant\n"; }
			print "$op: 0.0\n";
			next IN;
		}

		if ((abs((($data1[2] - $data2[2]) / &max($data1[2], $data2[2])) * 100)) < 10) {
			if ($debug) { print "latency is very close\n"; }
			print "$op: 0.0\n";
			next IN;
		}
	}

	# PHASE 2: group analysis
	$currdiff = 0.0;
	if ($phase2) {
		@res1 = &groupanal(@data1);
		@res2 = &groupanal(@data2);

		if ($res1[0] != $res2[0]) {
			if ($res1[0] > $res2[0]) {
				$more = \@res1; $less = \@res2;
			} else {
				$more = \@res2; $less = \@res1;
			}

			for($i = 1; $i <= $$more[0]; $i++) {
				if ($i <= $$less[0]) {

				} else  {
					if (($$more[$i]->{START} > $$less[$$less[0]]->{END}) && ($$more[$i]->{GROUPTOT} < 5)) {
						if ($debug) { print "ignoring small end group: $$more[$i]->{START} - $$more[$i]->{END}\n"; }
						$$more[0]--;
						for ($j = $i; $j < $$more[0]; $j++) {
							$$more[$j] = $$more[$j+1];
						}
					}
				}
			}
		}

		if ($res1[0] != $res2[0]) {
			if ($debug) {
				printf("$op: number of groups differ by %d\n", abs($res1[0] - $res2[0]));
			}
			print "$op: 100.0\n";
			next IN;
		} else {
			for($i = 1; $i <= $res1[0]; $i++) {
				$latdiff = abs((($res1[$i]->{AMOUNT} - $res2[$i]->{AMOUNT}) / &max($res1[$i]->{AMOUNT}, $res2[$i]->{AMOUNT})) * 100);
				if ((abs($res1[$i]->{START} - $res2[$i]->{START}) > 1) && ($latdiff > 50)) {
					if ($debug) {
						printf("$op: start of group %d differs by %d\n", $i, abs($res1[$i]->{START} - $res2[$i]->{START}));
					}
					print "$op: 100.0\n";
					next IN;
				} elsif ((abs($res1[$i]->{END} - $res2[$i]->{END}) > 1) && ($latdiff > 50)) {
					if ($debug) {
						printf("$op: end of group %d differs by %d\n", $i, abs($res1[$i]->{END} - $res2[$i]->{END}));
					}
					print "$op: 100.0\n";
					next IN;
				} elsif ((abs($res1[$i]->{END} - $res2[$i]->{END}) > 2) && ($i == $res1[0])) {
					if ($debug) {
						printf("$op: end of group %d differs by %d\n", $i, abs($res1[$i]->{END} - $res2[$i]->{END}));
					}
					print "$op: 100.0\n";
					next IN;
				} else {
					$currdiff += $latdiff;
				}
			}
		}
	}

	# PHASE 3: do a specific technique on the left-over profiles
	if (($command eq "grouplat") || ($command eq "groupops")) {
		$result = $currdiff;
	} elsif ($command eq "chisquare") {
		$result = (1 - &chstwo(\@data1, \@data2)) * 100;
	} elsif ($command eq "totlat") {
		$result = (abs($data1[2] - $data2[2]) / &max($data1[2], $data2[2])) * 100;
	} elsif ($command eq "totops") {
		$result = (abs($data1[1] - $data2[1]) / &max($data1[1], $data2[1])) * 100;
	} elsif ($command eq "correlation") {
		$result = &correlate(\@data1, \@data2);
	} elsif ($command eq "earthmover") {
		$result = &earthmover(\@data1, \@data2);
	} elsif ($command eq "quadratic") {
		$result = &quadratic(\@data1, \@data2);
	} elsif ($command eq "none") {
		$result = 50.0;
	} else {
		print "ERROR: unknown command $command\n";
		&usage();
		exit 1;
	}
	printf("$op: %f\n", &min(100.0, $result));
}

close INFILE1;
